转录本融合位点上下游序列获取

首先需要理解：基因，转录本(transcripts,isoform，mRNA序列)、EXON区域，cDNA序列、UTR区域，ORF序列、CDS序列 这些概念，一个基因可以转录为多个转录本，真核生物里面每个转录本通常是由一个或者多个EXON组成，能翻译为蛋白的EXON区域是CDS区域，不能翻译的那些EXON的开头和结尾是UTR区域，翻译区域合起来是ORF序列，而转录本逆转录就是cDNA序列。

STAR-Fusion 结果解读

通常建议大家对RNA-seq数据使用 STAR-Fusion 来检测转录本融合现象，得到的结果如下：

文件 'star-fusion.fusion_predictions.abridged.tsv', with the following format:

#FusionName           JunctionReadCount  SpanningFragCount  SpliceType           LeftGene                        LeftBreakpoint    RightGene                        RightBreakpoint   LargeAnchorSupport  FFPM        LeftBreakDinuc  LeftBreakEntropy  RightBreakDinuc  RightBreakEntropy  annots
THRA--AC090627.1      27                 93                 ONLY_REF_SPLICE      THRA^ENSG00000126351.8          chr17:38243106:+  AC090627.1^ENSG00000235300.3     chr17:46371709:+  YES_LDAS            23875.8456  GT              1.8892            AG               1.9656             ["CCLE","FA_CancerSupp","INTRACHROMOSOMAL[chr17:8.12Mb]"]

可以看到列比较多，值得我们关心的是转录本融合的左右两个转录本的ID及其融合位点在基因组的坐标，如下：

THRA^ENSG00000126351.8          chr17:38243106:+  
AC090627.1^ENSG00000235300.3     chr17:46371709:+  

上面的坐标可以无限多，文件命名为 pos.txt 吧，后面写脚本需要用到，值得注意的是作者软件举例应该是hg19的基因组坐标，目前几乎都是hg38了，所以后面我脚本也是基于hg38的。

这个时候我们可能需要设计引物来对该融合转录本 进行验证，所以会需要这个融合点左右两个基因的指定转录本的cDNA序列。

我们很容易拿到各个转录本的基因组坐标，但是融合点的基因组坐标不能简单对应到转录本cDNA序列里面坐标。我们的突破点，就是找到融合点的基因组坐标到底对应到转录本cDNA序列的哪个位置。

物种的基因及其注释文件的下载

这里首选：https://asia.ensembl.org/info/data/ftp/index.html

wget ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/cds/Homo_sapiens.GRCh38.cds.all.fa.gz
wget ftp://ftp.ensembl.org/pub/release-95/gtf/homo_sapiens/Homo_sapiens.GRCh38.95.chr.gtf.gz

首先解析基因组注释文件，就是gtf的：

 1       havana  exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1";

格式化代码如下：

zcat Homo_sapiens.GRCh38.95.chr.gtf.gz |grep -w havana |perl -F"\t" -alne  '{next if $F[2] ne "exon";@a=split(/\"/,$F[8]); print join("\t",$a[1],$a[5],$a[9],$F[0],$F[3],$F[4],$F[4]-$F[3]+1) }' > g2t2exon.txt
zcat Homo_sapiens.GRCh38.95.chr.gtf.gz |grep -w havana |perl -F"\t" -alne  '{next if $F[2] ne "CDS";@a=split(/\"/,$F[8]); print join("\t",$a[1],$a[5],$a[9],$F[0],$F[3],$F[4],$F[4]-$F[3]+1) }' > g2t2cds.txt

结果如下：

ENSG00000223972    ENST00000456328 1   1   11869   12227   359
ENSG00000223972    ENST00000456328 2   1   12613   12721   109
ENSG00000223972    ENST00000456328 3   1   13221   14409   1189
ENSG00000223972    ENST00000450305 1   1   12010   12057   48
ENSG00000223972    ENST00000450305 2   1   12179   12227   49
ENSG00000223972    ENST00000450305 3   1   12613   12697   85
ENSG00000223972    ENST00000450305 4   1   12975   13052   78
ENSG00000223972    ENST00000450305 5   1   13221   13374   154
ENSG00000223972    ENST00000450305 6   1   13453   13670   218
ENSG00000227232    ENST00000488147 1   1   29534   29570   37

可以看到，一个基因会有多个转录本，然后每个转录本有多个外显子。不同的转录本的外显子有重叠。值得注意的是很多gtf里面记录的基因，在cDNA序列文件里面是不存在的，不过这个不影响我们的任务。

我上面两个代码分别得到的是外显子和CDS的坐标，后来发现CDS才是正确的， 因为并不是所有的外显子序列都会出现在cDNA序列里面。

首先基因组坐标转为转录本坐标

接下来需要写脚本把我们转录本融合位点那个基因组坐标，转为其转录本的相对坐标，这个时候普通的shell脚本已经无能为力，需要python或者perl这样的编程语言啦，就是把我们的 pos.txt 文件和自己制作的 g2t2exon.txt 关联起来，所以需要使用关联数组这个东西。

当然，也可以使用大家最擅长的R语言展示咯。

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('input.txt')
head(a)
b=read.table('g2t2cds.txt')
colnames(b)=c('gene','transcript','exon','chr','start','end','exon_length')
head(b)
library(stringr)
re=NULL
tmp = apply(a,1,function(x){
  # x=a[1,]
  g=str_split(x[1],'[.]')[[1]][1]
  pos=as.numeric(str_split(x[2],':')[[1]][2])
  info=b[b[,1]==g,] ## one gene might has more than 1 transcripts.
  lapply(split(info,info[,2]), function(y){
    ## each transcript 
    apply(y,1,function(z){
      ## each exon.  
      if(z[5] <= pos & z[6] >= pos ){
        #print(c(z,pos))
        new = sum(as.numeric(y[y[,3] < as.numeric(z[3]),7])) + pos - as.numeric(z[5])+1 
        print(c(z,pos,new))
        re <<- rbind(re,c(z,pos,new))
      }
    }) ## each exon  
  }) ## each transcript 

}) ## each position 
colnames(re)=c('gene','transcript','exon','chr','start','end','exon_length','pos','new_pos')

代码非常复杂，需要一定编程水平才能理解。

1 ENSG00000121879 ENST00000643187 22 3 179234094 179235098 1005 179234680 3712
2 ENSG00000074800 ENST00000646370 13 1 8861007 8861429 423 8861330 1738
3 ENSG00000074800 ENST00000647408 12 1 8861002 8861429 428 8861330 1683
4 ENSG00000105976 ENST00000397752  1 7 116672392 116672577  186 116672464 73
5 ENSG00000105976 ENST00000456159 1 7 116672390 116672577 188 116672464 75
6 ENSG00000146648 ENST00000420316  7 7 55154011 55154152 142 55154029 1010
7 ENSG00000146648 ENST00000455089  6 7 55154011 55154152 142 55154029 888
8 ENSG00000077782 ENST00000526570  1 8 38427921 38430820 2900 38428753 833
9 ENSG00000037474 ENST00000502932 2 5 6603869 6604276 408 6604266 502

如上所述，融合位点的基因组坐标，成功转为了其对应的转录本序列坐标。

比如 ENST00000643187 这个转录本的坐标是

ENST00000643187.1 cds chromosome:GRCh38:3:179148574:179235098:1 gene:ENSG00000121879.4 gene_biotype:protein_coding

而融合位点在 179234680 ， 如果纯粹的使用它减去转录本起始坐标后是 86106 , 包含了大量的intron序列，所以需要找到其精准的外显子坐标。

比如这里的第22个外显子坐标是 3 179234094 179235098 ， 得到 586 的长度，再加上这个转录本前面的所有CDS的长度之和，最后是 3712 ， 就是该融合位点的转录本坐标啦。

然后根据转录本坐标及转录本序列获取

这个代码也不简单，需要读取我们下载好的cDNA序列文件：

rm(list = ls())
options(stringsAsFactors = F)
load(file = 're.Rdata')
fa=readLines('Homo_sapiens.GRCh38.cds.all.fa.gz')
fa=paste0( fa ,collapse = '')
fa=strsplit(fa,'>')[[1]]
fa=fa[-1]
tid=str_split(fa,'[.]',simplify = T)[,1]
seq=str_split(fa,']',simplify = T)[,2]
head(tid)
head(seq)
re=as.data.frame(re)
re=re[re$transcript %in% tid,]
re$tseq=seq[match(re$transcript,tid)]
head(re)
re$up=unlist(lapply(1:nrow(re),function(i){
  #i=1
  new_pos=as.numeric(re[i,9])
  l=nchar(re[i,10])
  if(new_pos>500){
    return(substring(re[i,10],new_pos-500,new_pos))
  }else{
    return(substring(re[i,10],0,new_pos))
  }

}))
re$down=unlist(lapply(1:nrow(re),function(i){
  #i=1
  new_pos=as.numeric(re[i,9])
  l=nchar(re[i,10])
  if((l-new_pos) >500){
    return(substring(re[i,10],new_pos,new_pos+500))
  }else{
    return(substring(re[i,10],new_pos,l))
  }

}))

判断前面转换好的转录本坐标是否扩充500后会越界，然后选取扩充500bp的序列即可。当然， 这个结果请一定要去IGV检查哦，我这里只是提供一种思路。

值得注意的是转录本其实还有正负链信息是需要考虑的。

STAR-Fusion 提供工具组建融合后的转录本

其实并不需要自行写脚本进行探索，研究一下 STAR-Fusion '--denovo_reconstruct' parameter. This requires that you include the '--FusionInspector inspect|validate' setting.

两个参数稍微有点区别

• If FusionInspector 'inspect' mode is invoked, then only the fusion-evidence reads are de novo assembled.

• If FusionInspector 'validate' mode is selected, then all reads aligned to the fusion gene contigs are assembled.

同理，上述代码可以使用perl或者python来完成， 小伙伴赶快试试看吧！

历史题目：

ps：历史题目的单个推文都可以点击阅读原文可以查看其他同学提交的作业！